home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / correlate.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  196 lines
  1. ;$Id: correlate.pro,v 1.8 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       CORRELATE
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the linear Pearson correlation coefficient
  11. ;       of two vectors or the Correlation Matrix of an M x N array. 
  12. ;       Alternatively, this function computes the covariance of two vectors 
  13. ;       or the Covariance Matrix of an M x N array.
  14. ;
  15. ; CATEGORY:
  16. ;       Statistics.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = Correlate(X [,Y])
  20. ;
  21. ; INPUTS:
  22. ;       X:    A vector or an M x N array of type integer, float or double.
  23. ; 
  24. ;       Y:    A vector of type integer, float or double. If X is an M x N 
  25. ;             array, this parameter should not be supplied.
  26. ;
  27. ; KEYWORD PARAMETERS:
  28. ;       COVARIANCE:    If set to a non-zero value, the sample covariance is 
  29. ;                      computed. 
  30. ;
  31. ;       DOUBLE:        If set to a non-zero value, computations are done in
  32. ;                      double precision arithmetic.
  33. ;       
  34. ; RESTRICTIONS:
  35. ;       If X is an M x N array, then Y should not be supplied; 
  36. ;       Result = Correlate(X)
  37. ;
  38. ; EXAMPLES:
  39. ;       Define the data vectors.
  40. ;         x = [65, 63, 67, 64, 68, 62, 70, 66, 68, 67, 69, 71]
  41. ;         y = [68, 66, 68, 65, 69, 66, 68, 65, 71, 67, 68, 70]
  42. ;
  43. ;       Compute the linear Pearson correlation coefficient of x and y.
  44. ;         result = correlate(x, y)
  45. ;       The result should be 0.702652
  46. ;
  47. ;       Compute the covariance of x and y.
  48. ;         result = correlate(x, y, /covariance)
  49. ;       The result should be 3.66667
  50. ;  
  51. ;       Define an array with x and y as its columns.
  52. ;         a = [transpose(x), transpose(y)]
  53. ;       Compute the correlation matrix.
  54. ;         result = correlate(a)
  55. ;       The result should be [[1.000000,  0.702652]
  56. ;                             [0.702652,  1.000000]]
  57. ;
  58. ;       Compute the covariance matrix.
  59. ;         result = correlate(a, /covariance)
  60. ;       The result should be [[7.69697,  3.66667]
  61. ;                             [3.66667,  3.53788]]
  62. ;
  63. ; PROCEDURE:
  64. ;       CORRELATE computes the linear Pearson correlation coefficient of
  65. ;       two vectors. If the vectors are of unequal length, the longer vector
  66. ;       is truncated to the length of the shorter vector. If X is an M x N 
  67. ;       array, M-columns by N-rows, the result will be an M x M array of 
  68. ;       linear Pearson correlation coefficients with the iTH column and jTH 
  69. ;       row element corresponding to the correlation of the iTH and jTH 
  70. ;       columns of the M x N array. The M x M array of linear Pearson 
  71. ;       correlation coefficients (also known as the Correlation Matrix) is 
  72. ;       always symmetric and contains 1s on the main diagonal. The Covariance
  73. ;       Matrix is also symmetric, but is not restricted to 1s on the main
  74. ;       diagonal. 
  75. ;
  76. ; REFERENCE:
  77. ;       APPLIED STATISTICS (third edition)
  78. ;       J. Neter, W. Wasserman, G.A. Whitmore
  79. ;       ISBN 0-205-10328-6
  80. ;
  81. ; MODIFICATION HISTORY:
  82. ;       Written by:  DMS, RSI, Sept 1983
  83. ;       Modified:    GGS, RSI, July 1994
  84. ;                    Added COVARIANCE keyword.
  85. ;                    Included support for matrix input.  
  86. ;                    New documentation header.
  87. ;       Modified:    GGS, RSI, April 1996
  88. ;                    Included support for scalar and unequal length vector
  89. ;                    inputs. Added checking for undefined correlations and
  90. ;                    support of IEEE NaN (Not a Number). 
  91. ;                    Added DOUBLE keyword.
  92. ;                    Modified keyword checking and use of double precision.
  93. ;-
  94. FUNCTION DATOTAL, Arg, Double = Double
  95.  
  96.   Type = SIZE(Arg)
  97.   if N_ELEMENTS(Double) eq 0 then $
  98.     Double = (Type[Type[0]+1] eq 5)
  99.   ArgSum = TOTAL(Arg, Double = Double)
  100.   if Type[Type[0]+1] eq 5 and Double eq 0 then $
  101.     RETURN, FLOAT(ArgSum) else RETURN, ArgSum
  102.  
  103. END
  104.  
  105. FUNCTION Cov_Mtrx, X, Double = Double
  106.   Nx = SIZE(x)
  107.   if N_ELEMENTS(Double) eq 0 then Double = (Nx[Nx[0]+1] eq 5) 
  108.   if Double eq 0 then one = 1.0 else one = 1.0d
  109.   if nx[0] eq 0 or nx[0] eq 1 then RETURN, one $
  110.   else begin 
  111.     VarXi = X - (X # REPLICATE(one/nx[2], nx[2])) # REPLICATE(one, nx[2])
  112.     VarXi = (VarXi # TRANSPOSE(VarXi)) / (nx[2]-1)
  113.     if Double eq 0 then RETURN, FLOAT(VarXi) else RETURN, VarXi
  114.   endelse
  115. end
  116.  
  117. FUNCTION CRR_MTRX, X, Double = Double
  118.   Nx = SIZE(x)
  119.   if N_ELEMENTS(Double) eq 0 then Double = (Nx[Nx[0]+1] eq 5)
  120.   if Double eq 0 then one = 1.0 else one = 1.0d
  121.   if Nx[0] eq 0 or Nx[0] eq 1 then RETURN, one $
  122.   else begin
  123.     VarXi = x - (x # REPLICATE(one/Nx[2], Nx[2])) # REPLICATE(one, Nx[2])
  124.     SS = VarXi^2 # REPLICATE(one, Nx[2])
  125.     SS = SS # SS
  126.     if Double eq 0 then Tol = 1.0e-6 else Tol = 1.0d-12
  127.     iSS = WHERE(ABS(SS) le Tol, nSS) ;Zero denominator signals undefined
  128.                                      ;Correlation.
  129.     if nSS ne 0 then begin
  130.       SS[iSS] = one
  131.       cm = (VarXi # TRANSPOSE(VarXi))/SQRT(SS)
  132.       cm[iSS] = !VALUES.F_NAN
  133.       if Double eq 0 then RETURN, FLOAT(cm) else RETURN, cm
  134.     endif else $
  135.       if Double eq 0 then RETURN, FLOAT((VarXi # TRANSPOSE(VarXi))/SQRT(SS)) $
  136.       else RETURN, (VarXi # TRANSPOSE(VarXi))/SQRT(SS)
  137.   endelse
  138. end
  139.  
  140. FUNCTION Correlate, X, Y, Covariance = Covariance, Double = Double
  141.  
  142.   ON_ERROR, 2  ;Return to caller if an error occurs.
  143.  
  144.   if N_PARAMS() eq 2 then begin  ;Vector inputs.
  145.  
  146.     Sx = SIZE(x)  &  Sy = SIZE(y)
  147.     if N_ELEMENTS(Double) eq 0 then $
  148.       Double = (Sx[Sx[0]+1] eq 5) or (Sy[Sy[0]+1] eq 5)
  149.  
  150.     Nx = Sx[Sx[0]+2]
  151.     Ny = Sy[Sy[0]+2]
  152.  
  153.     ;If not of equal length, take shortest.
  154.     if Nx lt Ny then begin
  155.       yMem = Y[Nx:*] & Y = Y[0:Nx-1]
  156.     endif else $
  157.     if Ny lt Nx then begin
  158.       xMem = X[Ny:*] & X = X[0:Ny-1]
  159.     endif
  160.  
  161.     ;Means.
  162.     sLength = MIN([Nx, Ny])
  163.     xMean = DATOTAL(X, Double = Double) / sLength
  164.     yMean = DATOTAL(Y, Double = Double) / sLength
  165.  
  166.     ;Deviations.
  167.     xDev = X - xMean
  168.     yDev = Y - yMean
  169.  
  170.     ;If not of equal length, restore input vector.
  171.     if Nx lt Ny then Y = [Y, yMem] else $
  172.     if Ny lt Nx then X = [X, xMem] 
  173.  
  174.     if KEYWORD_SET(Covariance) eq 0 then begin  ;Correlation.
  175.       Denom = DATOTAL(xDev^2, Double = Double) * $
  176.               DATOTAL(yDev^2, Double = Double)
  177.       if Double eq 0 then Tol = 1.0e-6 else Tol = 1.0d-12
  178.       iDenom = WHERE(ABS(Denom) le Tol, nDenom) ;Zero denominator signals 
  179.                                                 ;undefined Correlation.
  180.       if nDenom ne 0 then begin
  181.         Denom[iDenom] = 1.0
  182.         cr = DATOTAL(xDev * yDev, Double = Double)/SQRT(Denom)
  183.         cr[iDenom] = !VALUES.F_NAN
  184.         RETURN, cr
  185.       endif else RETURN, DATOTAL(xDev * yDev, Double = Double)/SQRT(Denom)
  186.     endif else begin ;Covariance.
  187.       if sLength eq 1 then RETURN, !VALUES.F_NAN $
  188.       else RETURN, DATOTAL(xDev * yDev, Double = Double) / (sLength-1)
  189.     endelse
  190.   endif else begin ;Array input.         
  191.     if KEYWORD_SET(Covariance) eq 0 then RETURN, $
  192.       CRR_MTRX(X, Double = Double) $  ;Correlation Matrix.
  193.     else RETURN, COV_MTRX(X, Double = Double)  ;Covariance Matrix.
  194.   endelse
  195. END
  196.